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In this work we develop an approach to treat correlated many-electron dynamics, 
dressed by the presence of a finite-temperature harmonic bath. Our theory combines 
a small polaron transformation with the second-order time-convolutionless master 
equation and includes both electronic and system-bath correlations on equal 
footing. Our theory is based on the ab-initio Hamiltonian, and thus well-defined 
apart from any phenomenological choice of basis states or electronic system- 
bath coupling model. The equation-of-motion for the density matrix we derive 
includes non-Markovian and non-perturbative bath effects and can be used to 
simulate environmentally broadened electronic spectra and dissipative dynamics, 
which are subjects of recent interest. The theory also goes beyond the adiabatic 
Born-Oppenheimer approximation, but with computational cost scaling like the 
Born-Oppenheimer approach. Example propagations with a developmental code 
are performed, demonstrating the treatment of electron-correlation in absorp- 
tion spectra, vibronic structure, and decay in an open system. An untransformed 
version of the theory is also presented to treat more general baths and larger systems. 
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I. INTRODUCTION 



The small-polaron transformation of the electronic Hamiltonian was originally developed 
in the 1960'sP^, and more recently revivecP^^ in the many-electron context. It is a classic 
example^ of the utility of canonical transformations in quantum physics. Its usefulness is 
well-established^ yet it is also experiencing renewed interestP^. In particular, second-order 
master equations in the polaron frame afford good results in all bath strength regimes^HlflU 
employing the variational technique of Harris and Silbey^ . In the many-electron case, 
there has also been some recent pioneering work toward developing random-phase approxi- 
mation equation d 92 * 93 l The electronic structure community has also produced some related 
wor kpI3HlQlD including phenomenologically damped-^ response theory. 

In this paper we formulate and implement a correlated many-electron master equation 
that overcomes several limitations of the adiabatic Born-Oppenheimer approximation, and 
includes effects such as excited-state lifetimes and vibronic structure. The basic goal of 
our method is to produce electronic spectra for small molecules that include the effects of 
coupling to an environment. With continued development our formalism could describe 
environmental localization of electronic states and the decay of quantum entanglement in 
correlated electronic systems. 

The theory we present exploits the polaron transformation to combine both electron corre- 
lation and system-bath couplings in a single perturbation theory. In the transformed frame, 
high-rank quantum expressions are dressed by environmental factors, which cause them to 
decay during dynamics. This introduces the possibility for environmentally induced decay 
of the correlations in an electronic system, thereby making the problem computationally 
more tractable. We also discuss how a general one-particle perturbation to the electronic 
system may be treated in a closely related untransformed version of the theory. Physically 
relevant coupling models of this form are numerous, and several examples include nuclear 
motion coupled to electronic degrees of freedom^ 110 * 111 !, Coulomb coupling to a nearby nano- 
particle surface^', the electromagnetic vacuumpS anc l perhaps even Coulomb coupling to 
surrounding molecules in a condensed phase. 

The present method is distinguished from previous work by a few characteristic features. 
Unlike virtually all master equation approaches, it treats the dynamics without assum- 
ing the many-body problem of the electronic system to be solved. We emphasize 
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that our correlated theory is one half of a final theory with the other involving a recipe to 
calculate system specific couplings with Ehrenfest type dynamics or another scheme. The 
reason being that, in general, the appropriate basis of a few electronic states to prepare a 
master equation can't be found a priori. As a result, this theory begins from a system-bath 
Hamiltonian which is well-defined and atomistic in terms of single-electron states and em- 
ploys a non-Markovian equation of motion (EOM) in place of phenomenological damping. 
Conveniently, we find that high-rank operator expressions responsible for the computational 
intractability of exact, closed, many-particle quantum mechanics are multiplied by factors 
which exponentially vanish in many circumstances. 

The picture of electronic dynamics offered by a master equation is complementary to the 
time-dependent wave-packet picture of absorptiori^HI2Q] j n latter approach the elec- 
tronic degrees of freedom are described as a superposition of a few adiabatic surfaces, and 
nuclear wave-packets are studied in the dynamics. In our approach the electronic degrees of 
freedom are described in all-electron detail, but the dynamics of the nuclei are approximated 
by the choice of the Holstein Hamiltonian. It is possible to combine a wave-packet calcula- 
tion of the bath correlation function with the untransformed equation of motion presented 
in this work, to leverage the physical strengths of both approaches. 

We demonstrate the implementation of the formalism in a pilot code, and apply it to 
calculate some dynamic properties of model small molecules. The spectra produced by the 
electron correlation theory are shown to compare favorably to related methods in the adia- 
batic limit. Vibronically progressed spectra are shown to be produced by the dressed theory, 
and a Markovian version is applied to a basic model of electronic energy transport between 
chromophores. Finally the undressed, uncorrelated theory is used to simulate the ultraviolet 
absorption spectrum of 1,1-diflouroethylene and compared with available experimental data. 



II. THEORY 
A. Hamiltonian 

As in the work of reference d 92 * 93 !, we use the non-relativistic ab-initio many-electron Hamil- 
tonian with a Holstein-type^ 2 ^ (linear) coupling to a bath of non-interacting bosons (summa- 
tion over repeated indices is implied throughout this paper unless stated otherwise), given 
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by 



H = F + V + H hoson + -ffei-boson = fpa\a p + V p °a ] s a\a p a q + u k b\b k + a\a v M v k (b\ + b k ) . (1) 

Here, a\ creates an electron in the single-particle basis state s, while b k creates a bosonic 
bath particle in the kth mode. F is the Fock operator of the reference determinant, and V 
is the two-electron part of the electronic Hamiltonian in the same single-particle basis. The 
third term is the boson Hamiltonian, Hf, and the last is the coupling of the electronic system 
to the bath. For a general bath mode with dimensionless displacement Q k , the bi-linear 
coupling constant M%, is related to the derivative of the orbital energy via M% = uo^ 
(no summation over p and q implied). Assuming this sort of coupling is only appropriate 
for nuclear configurations near the minimum of the Born-Oppenheimer well. There are 
prescriptions for defining these parameters from normal-mode analysis^ and molecular 
dynamics^' for the semiclassical treatment of anharmonicity. The latter approach can also 
include electrostatic fluctuations. For the purposes of this work we are only interested in 
the properties of the method we develop, and do not present any calculations of the bath 
Hamiltonian. In the following discussion, we will assume the one-electron parts of H and F 
to be diagonal in the (canonical) one-electron basis with eigenvalues e p = fg. 
We now introduce the displacement operator, 

exp[S] = exp[a p a p Ml{b\ - b k )\ (2) 
M k p = M k p /uj k , (3) 

which generates the polaron transformation. Since the operator S is anti-Hermitian, it 
generates a unitary transformation of the electron-boson Hamiltonian H — >■ H = e~ s He s . 
Explicitly, the polaron transformed Hamiltonian is given by 

H = F + H int} (4) 

where the one-particle part of the transformed Hamiltonian is given by 

F = F e ie + -ffboson = (c p - A p )aj,a p + u k b[b k , (5) 
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and we have introduced the reorganization energy X p = ^2 k M k /u>k- The transformed 
electron-electron interaction is given by 

H int = Vpq S alala p a q XlXlX p X g , (6) 

where the transformed matrix elements are 

V£ = V™ - 2(oo k M r k M s k )6 ps 6 qr (l - 5 sr ), (7) 

and we have defined the bath operators X p as 

X p = exp[M p k (bl-b k )}. (8) 

For future reference it will also be useful to define dressed electronic creation and annihilation 
operators, denoted by a\ = a\X\ and a s = a s X s . 

The key feature of the polaron transformation is that H has no electron-phonon coupling 
term. As a result, the two-electron and electron-boson parts of the original Hamiltonian 
are combined into a single term which now couples two dressed electrons and two dressed 
bosons (Eq. [6]). One should intuitively imagine the situation depicted in Fig. [TJ where two 
displaced electronic energy surfaces are dragged into alignment by the polaron transfor- 
mation, altering the coupling region between them, which now absorbs the electron-boson 
coupling. 

In what follows, our expressions will be derived in the interaction picture with respect to 
Eq. [6] and then switched to the Schrodinger picture. The harmonic nature of the bosons 
means that the correlation functions of X operators (in any combination and at multiple 
times) can be given as simple functions of u p , M and (5 = k^T^. The transformed electron- 
boson problem takes the form of the usual many-electron problem with a time-dependent 
electron-electron interaction. This allows us to harness powerful methodologies that stem 
from two distinct areas of research: quantum master equations^' and quantum chemistry 
methods^!. We exploit this feature to produce a model of electronic dynamics which treats 
system-bath dynamics and correlation effects within the same perturbation theory. 
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FIG. 1: A schematic representation of the polaron transformation. One-electron energies of 
the underlying Hamiltonian take the form of displaced parabolas coupled to the Coulomb 
interaction V. In the transformed frame, the boson effects are absorbed into a V, which 
depends on the boson operators, X. This makes V effectively time-dependent via the bath 

correlation function B(t). 

B. Electronic dynamics 

Having reviewed the polaron transformation in the previous section, we now combine 
it with the time-convolutionless perturbation theory (TCLj^^^ to arrive at the central 
theoretical result of this manuscript. Our goal is to derive an equation-of- motion for a dressed 
particle-hole excitation operator, which we denote b % a = o^a\ai (here i, j, k... are zeroth-order 
occupied levels and a,b,c... are unoccupied). We consider only the particle-hole part of the 
density matrix, commonly referred to as the Tamm-Dancoff approximatiorP^^TDA) , to 
simplify the derivation of our equations and avoid the possible issues associated with the 
non-linearity of the other blocks^'. To derive Fock-space expressions^', it is convenient 
to assume the initial equilibrium state is the canonical Hartree-Fock (HF) determinant i.e. 
the initial density matrix is | v I / (0))(\I/(0)| ~ |0)(0|, with |0) the HF determinant. This 
assumption relies on two approximations: 1) The first excited-state energy of the systems 
we will study is much larger than k^T, so that the system is effectively in the ground-state. 
2) The initial state of the system is weakly correlated and hence dominated by a single Slater 
determinant. Assumption 1 is clearly valid for the small molecular systems we will study, 
while assumption 2 requires more care and the effects of initial correlations are discussed in 
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detail in section II. E. We define projection operators: 

VA = V F (Tr b (A)) x (Xl..X) b , eg , and Q = 1 - V. (9) 

Vf is a Fock space projector onto maps between single-electron {a^a} operators like the 
typical projector of a partitioned electron-correlation perturbation theory ie: 

Vf{C\ ■ o q p a q a p — > r] s r a\a r ) = d (10) 
V F {C >=2 : o pq rs a\a\...a p a q ->• n/q^. a l> a l>--- a p' a q') = °- ( n ) 

This partitioning is consistent with the perturbative ordering of H by powers of V. We treat 
the description of the many-electron state in terms of only one-body operators (neglecting all 
higher density matrices) on the same footing as tracing over the bath degrees of freedomP^D, 
so both phenomena are easily incorporated in the same master equation. 

The effective Liouvillian becomes time-dependent due to the polaron transformation. 
Consequently, there exists no simple analytical formula for its Fourier transform. Instead, 
we must give a differential representation of the EOM for a particle-hole excitation, which 
can then be integrated numerically. Given these projectors, the time-convolutionless per- 
turbation theory^S (TCL) over the space of V isP^l. 

jVd{t) = (vC(t)V + dsVC(t)QC(s)V^J d{t) + X{t) (12) 

This is written in the interaction picture where £ = —i[V(t), ■]. The first term above 
is the uncorrelated part of the evolution!^, the second is a homogenous term reflect- 
ing correlation between the system and the bath, and the last term is an inhomogene- 
ity reflecting correlations of the initial state. The interaction picture perturbation^^ is 
V(t) = (\>V A -^ata r a s )x(Xt(t)Xt(t)X r (t)X s (t)),whereAj^ = (e 7 +e 6 +...-e a -ep-...). 
Expanding the first term over the one-particle space and moving into the Schrodinger picture 
we obtain (left arrows, are used to indicate the contribution of a term to dd/dt ): 

-^W(t) <- -i(e u -e v + V:?B; s u (t,t))5 s p (t) (13) 
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We have introduced a shorthand for the correlation function (lower time, upper 

time): 

B™° s P (t, s) = (Ai( S )JCl( a )X o ( S )X p ( fl )Aj(t)Aj(*)X P (0X.(0.> (14) 



Term (13) is comparable in dimension and physical content to the response matrices of 
configuration interaction singles (CIS)P^ with an attached, time-local boson correlation 
function. Because all arguments have the same time-index, B only applies a real factor 
(~ 1 as M — > 0) to values of the interaction and thus introduces no new time-dependence. 
The indices of the boson correlation function, and their order, are simply read-off the V 
integral they multiply. 

The second homogeneous term introduces bath correlation functions between boson op- 
erators occurring at different times (t and s) according to 

jVblit) <- ^ds{VC{t)QC{s)Vo{t)) (15) 

= -KM Q[ f B%g(t, s)V x °t(s)ds, ol(t))) (16) 

J t 

Moving the electronic part into the Schodinger picture and rearranging this becomes: 

= -IK*, Q\V£, d n m }}e^e^ f B*%(t, s)e^>ds (17) 

J t 

Expanding the commutators in Eq. [TTj applying Wick's theorem to remove the many vanish- 
ing terms, and enforcing the connectivity constraint, one obtains many topologically distinct 
terms. In our implementation this is done automatically before the execution of a simula- 
tion. The terms are easily related to terms which occur in the expansion of the second-order 
Fermion propagator (SOPPA)^Hl38] and diagonalization-based excited-state theories like 
CIS(D)P21. Since we employ the time-convolutionless perurbation theory, the oscillating e lAt 
factors are different than those which occur in Rayleigh-Schro dinger perturbation theories 
(in the energy-domain denominator ^3^;) and warrants further study. We some terms from 



a hole — > particle excitation which is obtained by multiplying both sides of equation (27) on 
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FIG. 2: The particle-hole — > particle-hole skeleton BrandoW^ diagrams in the 
second-order homogeneous term, the first two come from the VtV s o term in the expanded 
double-commutator and third from V t oV s . The fourth occurs in V t oV s , V s oV t and V t V s o. 
The boson correlation function depends on any lines emerging from the ellipses which 
represent virtual electron scattering events at times t, s. 

the left with (a^a„)t and applying Wick's theorem. 

~oKt) <- C^ l(AS)1 f B h *%{t, s)e« A >ds (18) 

J to 

This term has six indices overall, but can be factorized as follows: 

5t(t) <- I a Am where: P c {t) = V b a d j V»fe^ f B $*{t, s)e^>ds (19) 

The calculation of scales fifth order with the number of single electron states, and linearly 
with the number of bath modes (which go into the calculation of B(t)). This low scaling 
with number of bath modes is inherited from the approach of Silbey^l, Jang^ and Nazir^Sl. 
The algebraic version of the second term in Fig. [2] is: 

ol(t) <- V£V$%e«*%< f Bj$(t, s)e^ a >ds (20) 

J to 

This term is sixth order with the size of the system, with the boson correlation function 
preventing a desirable factorization of the V-V contractions. 14 sixth-or-less-order terms are 
found which couple hole-particle excitations to each other; skeletons^ of these are shown 
diagrammatically in Fig. ^ with explicit expressions for each of these terms given in the 
supplementary informational. We should note that our formalism lacks several terms which 
occur in the SOPPA because of the Q projector and the absence of V t o vac oV s ordered terms 
that are invoked in the formal expansion of the time-ordered exponential. 

As described so-far the EOM is not time-reversible even when M = because neither 
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the homogenous term nor the normal excitation operator take the form of an anti-Hermitian 
matrix.'^!. We remedy this by adding the terms which result from the normal excitation's 
operator's Hermitian conjugate, ie: the indices of the vacuum and o are swapped, and the 
signs of A s and prefactor are flipped and the new term is added to the other perturbative 



terms. For example Term 20 becomes the following two terms: 



3£ <- ^^|^(M)(e ,(A - )(t - 8) )^ 



T ~. 



5- <- T V£V$%J B$£(t,8)(e« A #W)d8 (21) 

Making this modification, the linear response spectrum of the adiabatic model is not signif- 
icantly altered, but the adiabatic norm conservation is enforced, and changes by less than 1 
percent after 1700 atomic units in the case of if 4 with 4 th order Runge-Kutta and a timestep 
of 0.05 au. 



C. Dipole Correlation function 

Like any canonical transformation^^, the polaron transformation preserves the spectrum 
of the overall electron-phonon Hamiltonian but the statistical meaning of the state related to 
a particular eigenvalue is changed. In other words: the operators of our theory are different 
objects from the adiabatic Fock space. To lowest order in system-bath coupling^', electronic 
observables like the time-dependent dipole correlation function, which predicts the results 
linear optical experiments, can be obtained by generating the transformed property operator 
p, = e s fie~ s = ^a\ajXjXj. Within the Condon approximation the dipole correlation 
function is then the product of the electronic trace and bath trace over the X operators 
introduced in fi. This adds a bath correlation function to the usual observable. Taking all 
the relevant expectation values gives the following explicit numerical formula for the dipole 
moment expression: 

C d - d (t) = J2{ViaOia(t)n jb o jb (0)} ■ Tr B {Xl(t)X t (t)Xl(0)X 3 (0))} (22) 

ijab 
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This expression for the dipole-dipole correlation function assumes the bath remains at equi- 
librium. CIS likewise only offers a zeroth-order oscillator strengthens. The spectra in this 
work are generated by "kicking" the electronic system with the dipole operator instanta- 
neously, and Fourier transforming the resulting dipole-dipole oscillations. 

D. Untransformed Version 

Because of the polaron transformation the above formalism is accurate in the strong bath 
regime with diagonal system-bath couplings. To treat an off-diagonal, weak coupling one can 
develop the complementary untransformed theory. With a bi-linear system-bath coupling 
of the form: 

H sh ^ a \ a ^{bl + b k ), (23) 

development of an uncorrelated particle-hole equation of motion follow similarly to the 
above with H s b taking the place of V. The untransformed version also makes it possible to 
introduce an Ehrenfest scheme for the nuclear bath which may be pursued in future work. 
The projector of the untransformed version is simply the equilibrium trace over b operators. 
The TCL produces second-order contribution of the form: 

6{t)<r- [ [H sb (t),[H sb (s),d(t)]]ds (24) 

which after translation into the Schrodinger picture and application of Wick's theorem pro- 
duces several terms^SI similar to the following: 

6?(t) <- Mj: 1 o a m {t)e^ f e iA T°C(t - s)ds (25) 

J t 
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with t = t — s. The correlation function C m (t) = (6j n (t)6 m (0)} is given by the usual 
formula^! 



These contributions to the equations of motion scale with the 4 order of the system size; 
two-orders of magnitude cheaper than the transformed version. They can be combined with 
just the uncorrelated part of the particle-hole equation of motion, or added to the correla- 
tion terms developed above with the bath factor removed. To treat a continuous bath of 
oscillators, one can introduce a continuous parameterization of (M&) 2 = J(u)/u 2 = rji-^'^ 
as an ohmic spectral density under the assumption that the coupling to the continuous bath 
is the same for every state up to a factor. 

Furthermore, we find that the reorganization energy is given by Aj = uiid^ which can be 
subtracted from the original Hamiltonian. 

E. Initial Correlations 

In the above we have expanded on existing master equations by giving the electronic 
system a many-body Hamiltonian with Fermionic statistics. We have also improved typi- 
cal electronic response theories by incorporating bath dynamics, but we have assumed that 
the initial state was a single determinant and employed the corresponding Wick's theorem. 
These approximations are made in many other treatments of electronic spectra^ 10 * 148 * 149 !, 
and are acceptable for situations where a molecule begins in a nearly determinantal state. 
Optical absorption experiments of gapped small molecules belong to this regime. To rig- 
orously study electron transport in a biased junctionpSl, or otherwise more exotic initial 
state requires a treatment of initial correlations^^^, as in the non-equilibrium Green's 
function(NEGF) methocP^^. These methods first propagate an initial determinant in 
imaginary time^' thermalizing and correlating the system before the dynamics, but numer- 
ical applications of the NEGF formalism require storage an manipulation of a state variable 
with three indices: G pq {u), and are usually limited to small systems^!, and short times. We 





(26) 
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can similarly perform an imaginary time propagation, and have performed some exploratory 
calculations which are described in the supporting supporting material^. 

Instead of performing an imaginary time integration, one can replace the usual Wick's 
theorem, and build the theory beginning from a state which is already correlated. Extended 
normal ordering^!, makes it possible to take expectation values with a multi-configurational 
reference function via a generalization of Wick's theorem^. Using the extended normal or- 
dering technique, the expressions of the present paper can be promoted to treat a general 
initial state directly, so long as the density operator of that state is known. This approach 
would eliminate the need for any adiabatic preparation of the initial state, but would gen- 
erate more complex equations, and is a matter for future work. 

For weakly correlated systems it's reasonable to keep the approximation that the usual 
Wick's theorem applies to the initial state of the electronic system, and instead of adiabat- 
ically preparing an initial state, choose a perturbative approximation to the correlated part 
of the initial state, Qo(0) = J° oo dsC(s)o(0) 1 and treat the first inhomogenous^' term of 
the master equation: 

T(t)=VC(t)Q f dsC(s)o(0) (27) 

J-oo 

The bath parts of the inhomogenous term are known to be relatively unimportant from 
studies of tight-binding Hamiltonians^, and only slightly perturb the results of a propagation 
for short times, and so we will not include them. 

It is interesting to note that with the addition of this inhomogenous term, there is a 
correspondence between the present theory and a model of electronic linear response derived 
from an effective Hamiltonian, CIS(D)P^. The CIS(D) excited states are the eigenvectors of 



a frequency dependent matrix, A^f^- (u), with the dimension of the particle-hole space: 

A C JS D) H-H^-^0^ + VVQf 2 V (28) 

In the above T2 is the second-order excitation operator from Moller-Plesset perturbation the- 
ory. These terms correspond respectively to the Fourier transform of our VVV, VVQVV, and 
I(t) terms with different denominators. This correspondence suggests the present method 
should produce linear response spectra of quality similar to CIS(D), which is usually slightly 
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Approximation 



Extension 



TDA 

Orbital relaxation 



Second-Order in V 



(ji{t)) « Tr(fL ■ o(t)) 



^e q ) ~ |0) 



Extended normal ordering! 

TCL-4S3 

Use basis commuting with /i 
Include the additional blocks^S. 

Variational conditional 



TABLE I: Limitations of this work and how they may be relaxed. Orbital relaxation isn't 
really an approximation, per-se, but would be beneficial for the results of the perturbation 

theory. 

better than TDDFT. Because the formalism is somewhat involved and many approxima- 
tions have been made, we have summarized the limitations of this work in a table ([!]) with 
references that point to possible improvements. The formalism is now developed to the 
point where particle-hole excitations can be usefully propagated, and the main features of 
the approach can be demonstrated in calculations. 



To incorporate both bath and electron correlation effects it was necessary to write down 
a second-order, time-local EOM for electronic dynamics based on the time-convolutionless 
perturbation theory which we will call "2-TCL". The zeroth order poles of the correlation 
terms in this theory differ from those which occur in other second-order theories of electronic 



response (SOPPA^, ADC(2)P3S1 ; CIS(D)P2I, and CCSPH) which arise from perturbative 



partitioning of what is essentially an energy-domain propagator matrix. Interestingly, the 
denominator of the present theory is naturally factorized and in the adiabatic limit all terms 
can be evaluated in fifth order time, unlike Rayleigh-Scho dinger perturbation theory which 
requires a denominator factorization approximation to avoid a 6-index denominator^. To 
verify that the electronic part of this work is indeed a reasonable model of electronic dynam- 
ics and check our implementation (signs, factors etc.) it is useful to compare an adiabatic 
spectrum (a calculation with no bath coupling, M — Y 0) to one arising from exact diago- 
nalization. We've coded the above formalism into a standalone extension of the Q-Chem 
packagd^ from which we take the results of some other standard models. The particle hole 



III. RESULTS 



A. Adiabatic (M — > 0) spectrum 
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FIG. 3: Left: Adiabatic dipole absorption spectra of H 4 in the minimal basis, CIS which is 



a part of the first term in 27 CIS(D), and the present theory. Numbers indicate maxima 
in eV and horizontal axis is given in Hartree. The largest maxima have been normalized to 

the same value. Right: adiabatic 2-TCL applied to the BH 3 molecule. The placement of 
red peaks nearer to blue than green indicates success as a perturbation theory. The atomic 

geometries (given in a footnote, are also shown). 



equation of motion is integrated with the Runge-Kutta 4 th /§ th method with an adaptive 
time step. Propagations in all three directions are followed at once, and the resulting time 
dependent dipole tensor is Fourier transformed to produce the spectra presented below. 
Bath integrals are calculated with 3rd-order Gaussian quadrature with the same time-step 
as the electronic Runge-Kutta integration, or integrated analytically in the adiabatic limit. 
The exact results and moments shown below come from the PSI3^ program package. 

To check the adiabatic theory, we present calculations of dipole spectra on the H4 and 
BH 3 molecules^'. In both cases the molecules have been stretched from their equilibrium 
bond lengths to a geometry where correlation effects are stronger^' and excitations are 
anomalously low-in-energy because of near degeneracy of the single-particle levels. We 
prepared the minimal basis molecule of H4 in a density excited by the dipole operator, and 
propagated it for 250au. The dipole-dipole correlation functions C a p(t) = (fj, a (t)f/,p(0)) were 
collected during the simulation and Fourier transformed. The real part of this spectrum 
of spherically-averaged dipole oscillations is compared against stick spectra with the height 
given by the transition moment at the poles of an exact adiabatic calculation within this 
basis, and related theories. One sees that the propagator of the present work is a meaningful 
correction to CIS, moving poles away from the green positions towards the blue (Fig. [3]) as 
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CIS involves exact diagonalization of the single-excitation part of the many electron Hamil- 
tonian while the present theory involves correlation effects as well. In the adiabatic limit of 
the response model, only correlation effects (and not bath effects) are present. The excited 
state near 12.52eV in H4 has increased in energy from the green position of 11.54eV toward 
the exact pole at 13eV. A more strongly correlated state just below this is unresolved. This 
is likely because the dipole moment operator doesn't couple the reference determinant to 
this multi-configurational state. The higher energy peak of the spectrum is brought into 
good agreement with the exact result. Similar performance is seen in the case of a distorted 
BH 3 molecule, albeit with corrections of the higher energy peaks near 23eV being too small. 

B. Vibronic features 

The spectra of Markovian system-bath perturbation theories take the form of Lorentzians^^ 
at the poles of the response matrix. Because the present theory is non-Markovian it should 
be capable of yielding new poles off main transitions. To evaluate this effect we add a strong 
bath oscillator at 1600cm -1 and calculate absorption spectra. We choose this bath oscil- 
lator at high energy to observe vibronic peaks in a reasonable amount of propagation time 
(1700 a.u.), and choose a correspondingly high temperature (4397.25 K) so that there are 
several peaks in the progression which should take-on a Boltzmann-ian shape. The resulting 
absorption spectrum is pictured in Figure |4j with a close up of the promised progression. It 
is relevant to wonder whether the vibronic peaks are simply the result of the bath displace- 
ment operators present in the dipole correlation function expression, or the result of the 
polaron density matrix dynamics. Simply eliminating the bath-correlation function from 
the dipole expression and generating the same spectrum, a vibronic progression still results, 
indicating that it is the non-Markovian dynamics of the system operators which provides 
the vibronic progression. 

C. Energy Transport and Markovian Evolution 

The continuous integration of rank-6 bath correlation tensors required for the Non- 
Markovian propagation executed above is a rather costly proposition for large systems. This 
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FIG. 4: Polaron-transformed absorption spectrum of H4 with a Lorentizan bath oscillator 
at 1600cm -1 . Each adiabatic peak is split into a Boltzman-weighted vibronic progression 

magnified at inset. 



is especially true if one would like to study incoherent electronic energy transport which takes 
place in times on the order of picoseconds, roughly a million times more than the electronic 



timestep required to integrate Eq. (27) for a typical molecule. A useful approximation to 
overcome this is a Markov approximation, by which we mean taking the limit of the integral 
to infinity for each term, such as: 

K% = , lim f ! &- s)(e^^)ds (29) 

where R is now a factor replacing the integral expression. Each term possesses its own R 
tensor, but these must only be calculated once. Numerically, this limit can be taken in 
the case of a continuous super-ohmic bath with cutoff frequency u c by introducing a cutoff 
time t c , above which the bath correlation function is assumed to be equal to it's equilibrium 
value, which is a good approximation for (3u c (t c ) 2 » 2. 

We have used the letter R to suggest the analogy between this time-independent rate 
tensor and that occurring in the Redfield theory^USS although in this theory there are 28 
such fifth and sixth rank tensors. These tensors also differ because they use the transformed, 
rather than the bare correlation function. The value of the integral above only depends on 
the values of a Laplace-transformed B{uj) = J °° e lult B(t, 0) at the zeroth order electronic 
frequencies (A). The Markovian rate matrix can then be calculated once in sixth order 
time, giving effective kinetic rates which include the effects of correlation and bath coupling. 
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FIG. 5: Adiabatic linear response spectrum of a minimal model of transport between 
electronically excited states, two H2 molecules in the DZp^ basis. 



These can be used to calculate dynamics at a drastically reduced fourth order cost (since the 
perturbative EOM then takes the form of a single matrix product with time-independent 
effective Hamiltonian) or diagonalized directly to obtain spectra. The frequency independent 
nature of this term means that no new peaks can appear due to correlation or system- 
bath coupling. Only damping of dynamics between correlation and bath shifted CIS-like 
poles leading to Lorentzian spectra can be expected within the confines of this Markovian 
approximation. 

To give an example of how this methodology could be applied to transfer of electronic 
energy, we examine two Hydrogen molecules separated several bond lengths^' in a DZ 
basis. The basis has been expanded to allow for the induced dipole moment between the 
two molecules. Since these two molecules are well separated and held at nearly their Born- 
Oppenheimer minima there is no longer any strong correlation or spin-frustration in this 
system. The adiabatic linear response spectrum is shown in Fig ([5]), again showing good 
agreement with the stick spectra emerging from exact diagonalization as in if 4 . The CIS 
states of this pair of molecules are well-localized. The splitting between the two peaks around 
18eV corresponds to a weak coupling between them, and a distortion to the molecule on 
the right which has been imposed to lower the energy of the state localized on that side 
so that a bath can induce transfer of population. We initialize the pair of molecules into 
an even superposition of their excited states and couple a 2831cm -1 oscillator at 273. OK to 
the HOMO, HOMO-1, LUMO and LUMO+1 with dimensionless M's of (0.025, 0.05, .165, 
and .055) respectively. Over the course of ~60fs (Fig. pn the overlap of the time-dependent 
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FIG. 6: Transfer of probability between H2 molecules in the DZ basis. After being 
initialized in a transition density which is the superposition of the adiabatic states at 17.3 

and 18.1 eV, probabilities relaxation proceeds because of interaction with the bath. To 
zeroth-order in bath coupling State 1 is localized on the left, undistorted H 2 , and State is 

localized on the right, distorted H 2 . 

state with the higher energy state has halved the overlap with the lower energy state has 
grown, and the overall norm of the state has decreased by about 10 percent, corresponding 
to non-radiative decay. The rapid beating between the states corresponds to what would 
tight binding- model would describe as the coherence between collective states and 1. 
The evolution of this "coherence" isn't captured well in the Tamm-Dancoff approximation 
because we have truncated the blocks of the propagator coupling particle-hole excitations 
and their conjugate hole-particle modes. 

D. Untransformed and uncorrelated version 

The transformed version of the theory has several interesting formal advantages, but 
to treat larger systems or off-diagonal system-bath couplings we have also presented the 
untransformed version of our theory which is similar to a propagation of the CIS wave- 
function with dissipative bath terms. Even in our relatively rudimentary code it is possible 
treat much larger systems with the undressed version especially if the electron correlation 
terms are neglected. To demonstrate the idea we have simulated the ultraviolet photoabsorp- 
tion spectrum of 1,1-diflouroethylene^-^. We will compare the result of an untransformed 
calculation (Fig. [7| with available experimental dataP^ for the valence 71 — > 71* transition. 
Bath bosons are positioned at the experimentally raman active frequencies, the C=C stretch 
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FIG. 7: Experimental and simulated UV absorption cross-section of 1,1-diflouroethlyene in 
the region dominated by the it — > n* transition produced by 1500a.u. of propagation at 
303K. The untransformed version of the theory is employed and correlation terms are 

neglected. 



at 0.214eV, the CF 2 symmetric mode at .114eV, and the CH 2 rocking mode at .162eV. We 
have also added a super-ohmic spectral (n=3) density with coupling constant a = 0.01 and 
cutoff frequency of 5580 cm -1 . Because of a mixture of basis limitations and the absence of 
electron correlation the simulated absorption is centered an electron-volt higher in energy 
than the gas-phase experiment, but both transitions appear in a qualitatively similar way. 
The somewhat irregular vibronic progression of the experimental peak hints at limitations 
of approximating the nuclear degrees of freedom as a harmonic bath. This suggests that in 
future work it may be interesting to explore more accurate ways to calculate (6 t (t)6(0)), for 
example by propagating frozen gaussians or pursuing an Ehrenfest-type 1?:i scheme. In this 
untransformed version semiclassical or wave-packetP^ approaches to calculate (6^(t)6(0)) are 
easily combined into the electronic equation of motion. 

IV. CONCLUSION 

In this paper we have presented a model of correlated electronic dynamics that is trans- 
formed by a Bosonic bath. This allows us to study the effects of environmental dephasing 
on electronic excited states without assuming a reduced model for the electrons and their 
coupling to each other. Intriguingly, factors appear in the theory which damp high-rank 
operator expectation values. The possibility of exploiting this damping to reduce the 
cost of high-level electronic structure is interesting^!. The adiabatic limit of the prop- 
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agation developed here offers a useful correction to the placement of electronic energies. 
The linear-response spectrum of the complete model has been shown to exhibit vibronic 
structure. A Markovian, Tamm-Dancoff approximation to an energy transport problem has 
been examined, leading to bath-induced transfer of electronic energy between H2 molecules, 
although coherence decay cannot yet be captured because several off-diagonal blocks of 
the transition density matrix have been neglected in the EOM of this initial work. The 
untransformed version of this formalism could be made efficient enough to study non-linear 
optical experiments^H^ via propagation. 

Accurate protocols for providing a per-orbital spectral density must be developed and 
the accuracy of the resulting lineshapes and lifetimes must be assessed. Simply projecting 
the orbital energy gradient onto the normal modes of the molecule is an obvious choice, 
but it would be appealing to treat the influence of surrounding molecules in a similar way. 
The Tamm-Dancoff equation of motion provided in this work propagates only one-block of 
a much larger transition density matrix which must be treated to capture coherence phe- 
nomena between particles amongst themselves (and holes). The treatment of a thermalized 
initial condition is an interesting, but computationally demanding question. The payoff for 
developing these additional features would be electronic excited states with more realism 
than can be offered in the adiabatic picture. These states would be naturally localized, with 
a size depending on vibronic structure and temperature and evolve and relax irreversibly. 



V. APPENDIX: EXPRESSIONS FOR THE CORRELATION FUNCTIONS 
AND FACTORIZATION 

Time-ordered harmonic correlation functions (HCF's) of all orders of the type (XXt..XX^ 
were made available in Dahnovsky's pioneering workP^ by expanding the exponentials in X 
as power series in Ak = (bke ia>kt — b\e~ lu)kt ) and applying Wick's theorem to the resulting A 
operator strings paying careful attention to the combinatorial statistics. Since we employ 
a master equation theory rather than a Green's function theory, our version of the HCF 
depends on the sign of (t — s) but is otherwise the same after making the simplifications 
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which appear in our second order theory. 



B ££E£(M) = (xt i (t)xt 2 (t)x P3 (t)x P4 (t)xt i ( s )xt 2 ( s )x f/3 ( s )x g4 ( s )) = 

exp{- £ ^Coth(^ m /2)(M^ 3 ^(m)) 2 } • exp{- Y^^llt^)M^l{m))F m {t - s)} 



in 



where: F m (t) = Coth(/3w m /2)Cos(u; m t) — iSin(u m t) 

(30) 

We use an abbreviated notation (the same as A), M^--(s) = (M^ + M^...-M l s -M J s ...). The 
correlation function above includes equilibrium value of the HCF. The qualitative behavior 



of the real part of bath integrals in terms like 20 governs relaxation and is worth comment 



If u s ~ A > at time s, and M at t M at s < then summed over all terms Re(B(t)) is 
negative (causing relaxation) for a time on the order of l/(u s — A s ), after which it oscillates. 
If A < then relaxation occurs if M at t M at s > 0. In most applications of Master equations 
a single, positive, spectral density is assumed for all states which basically parameterizes 
M a t tM at s as a function of u. This approximation is more than a mere convenience; if M 
is assigned generally and different states are allowed to couple to a single frequency with 
different strengths it's quite easy to for some elements of the EOM to have M at tM at s > 
causing exponential growth in the Markovian limit and at short times. 

To treat a continuous number of bath oscillators one can introduce a continuous param- 
eterization of M called the spectral density, Ji(u), M^ a = J °° sj Ji{uS) j 'uS(u — uo a )duj. The 
renormalization of the electronic integrals is clear given this form. With a relatively simple 
functional form for J, such as a super ohmic spectral density with cutoff parameter u c , 
Jiiuj) = ^6^ e_W//aJc ! the time dependence of of correlation function B{t } s) can be analyti- 
cally calculated^ to a good approximation. So long as oo c is the same for all single-electron 
states, and only 7/j changes the whole formalism works identically with ^frjl taking the role of 
■ We note that the requirement that the correlation function be easily integrable is only 
a pre-requisite for the polaron transformation. If only the time-convolutionless equation of 
motion is used, any correlation function which is known can be easily incorporated. 

Without further approximation a sixth-order number of B's must be calculated and in- 
tegrated (in our code a third order Gaussian Quadrature with the electronic time-step was 
found sufficient), since at least two indices are shared between each V. Consider the cost 
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limiting term Eq. (20). It would be advantageous to evaluate the implied sum over c and 
make a 5-index intermediate, since then the remaining indices are also a fifth order loop, 
unfortunately the HCF depends on index c. On the other hand, High-rank operator strings 
are exponentially damped by the presence of this HCF in the Markovian limit. It seems 
likely that this feature could be used as a new locality principle which would lift the curse 
of dimensionality in the strong bath regime. 
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